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We report molecular dynamics calculations on two-component, cold (1.2 > T > 0.56 in natural 
units), three-dimensional Lennard-Jones fluids. Our system was small (7813 A, 7812 B particles), 
dense (N/V = 1.30), and distinctly finite (L x L x L cube, periodic boundary conditions, with 
L = 22.96(taa), o"aa being the range of the AA interaction in the Lennard-Jones potential J7y = 

4e ^(-^-) 12 — (~r") 6 j • We calculated spherical harmonic components Qlm for the density of particles 

in the first coordination shell of each particle, as well as their spherical invariants {Qt), time- 
correlation functions and wavelet density decompositions. The spherical invariants show that non- 
crystalline septahedral ((Qr)) ordering is important, especially at low temperature. While (Qio) 
could arise from icosahedral ordering, its behavior so closely tracks that of the nonicosahedral (Qii) 
that alternative origins for (Qlo) need to be considered. Time correlation functions of spherical 
harmonic components are bimodal, with a faster temperature-independent mode and a slow, strongly 
temperature-dependent mode. Microviscosities inferred from mean-square particle displacements are 
exponential in static amplitude (Q?), and track closely in temperature dependence the orientation 
density slow mode lifetime. Volume wavelet decompositions show that when T is reduced, the 
correlation length of (Q7) increases, especially below T = 0.7, but the correlation length of (Q|) is 
independent of T. 



I. INTRODUCTION 

Many fluid systems, on cooling below their freezing 
points, decline to crystallize, at least on experimentally 
accessible time scales. The fluids instead become glasses, 
highly viscous liquids or amorphous solids that appear 
rigid, but that lack long range crystalline order. De- 
spite enormous experimental, theoretical, and simula- 
tional study, the nature of the glass transition remains 
unclear. Is the transition toward the glass a dynamic 
effect, or does it reflect thermodynamic issues? A vast 
range of possible physical quantities might be correlated 
with glass formation, only a few of which have been stud- 
ied systematically. 

This paper treats orientational order, and fluctuations 
in orientational ordering, in cold Lennard-Jones fluids. 
Our results are based on molecular dynamics simulations. 
Orientational ordering in glass-forming liquids and other 
systems has been studied extensively pi H,H,H, OJ LJ0 H> 
The work reported here differs from earlier studies 
in that we (i) compute orientation correlations that cor- 
respond to non-crystalline ordering, (2) perform wavelet 
decompositions of the orientation correlation density, and 
(3) determine the temporal evolution of fluctuations in 
orientational ordering. 

There are a variety of ways to characterize angular 
order around a given particle. For example, one can 
measure the 'bond angle distribution', the relative likeli- 
hood that the displacement vectors from the given atom 



to two neighboring atoms have a particular value0. 
The bond angle distribution, which reflects an aspect 
of the three-particle distribution function g( 3 \ri, r2, r^), 
has recently been shown to be accessible to experimen- 
tal measurement 8J. Steinhardt, et al. 1] used a spher- 
ical harmonic decomposition of the near neighbors of a 
given particle, an idea referred back to the earlier work 
of Frank (11). Spherical harmonics are a complete set of 
functions over the sphere, so the angular distribution of 
(for example) an atom's nearest neighbors can always be 
expanded in spherical harmonics of their density. 

To characterize an angular distribution, one chooses an 
atom of interest, and a set of other atoms whose angular 
distribution around the atom of interest is to be charac- 
terized. The unit vector f j from the atom of interest to 
each atom i of the N other atoms are generated. The 
LM th spherical harmonic component Qlm of the f j is 
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where Ylm is the LM th spherical harmonic, and in the 
spherical harmonic fj is represented by the polar angles 
(6, 4>). The polar angles refer to a choice of the angular 
origins. To eliminate the dependence of Qlm on this 
nonphysical feature, one forms spherical invariants 
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that are independent of coordinate axesp]. 

Many prior treatments of particle orientations focused 
on a search for icosahedral, cubic, or other crystalline lo- 
cal ordering in the fluid, including local orderings that are 
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not crystallographic because they do not lead to space- 
filling packings. Partial crystalline orderings are charac- 
terized by particular values for some of the Qf For ex- 
ample, icosahedral order around a center leads to nonzero 
Qf for £ = 6, 10, 12 while face-centered-cubic and body- 
centered-cubic clusters have Qf ^ for £ = 4, 6, 8. As 
applied to liquids, the ensemble-average {Qf) determined 
from simulation or elsewise are compared with the {Qf) 
that would arise from ideal local order of one sort or an- 
other. For historical reasons, prior work has largely been 
limited to {Qf) for even £. However, the Qf are positive 
scmidcfinitc, and fluctuations can certainly lead to non- 
zero values for any Qi, so therefore it is inescapably the 
case that {Qf) is nonzero for every £. 

The spherical invariant Qf is determined by the instan- 
taneous particle positions, and is therefore a function of 
time and space: Qf changes as particles that are ini- 
tially nearby move or are replaced with other particles. 
At a given moment, Qf also varies from atom to atom. 
Wang and Stroud [T(| report studies of a spatial corre- 
lation function (QLM^a^LM^a + r )}for the spherical 
harmonic components around two atoms separated by 
r. The time correlation function of the time-dependent 
Qlm(*), namely 

4 1 

C?\t) = (Qm(0)Q* m {t)) (3) 

M=-l 

describes the relaxation of orientational order fluctua- 
tions around a given particle. Sanyal and Sood 4] stud- 

(2) 

ied a variation on C\ (t) , in which an elaborate weighted 
sum of Qlm over all particles in small cubic subsections 
of the system was taken. The time correlation function of 
the weighted sum was computed for each cubelet. Sanyal 
and Sood's 'bond-orientation correlation function' had a 
relaxation that was fit adequately well by a stretched ex- 
ponential in time. 

We further treat wavelet decompositions of the density 
of the {Qf). Wavelets [T^ provide via functional trans- 
formation a representation of a function g(x). A wavelet 
transform g(a, b) of a function g{x) is obtained as 

g(a,b) = J dx g(x)F(a,b;x) (4) 

where a and b are wavelet parameters, the dilation and 
the translation, respectively, and where the F(a,b;x) are 
the basis vectors of the wavelet transformation. The 
g{a,b) are the wavelet components of the function g(x). 
The basis vectors all follow from a single mother wavelet 
function f(x) via 

F(a,b;x) = f{ax-b); (5) 

the effectiveness of the transform being determined by 
the choices of /, a, and b. One choice of dilation 
and transformation variables relies on binary decimation, 
namely a = 2 m and b = n£o2 m , £ a being a basic length 



and n and m being integers. The g{a,b) allow one to 
regenerate the original g(x), namely 

g(x) = J da db g(a,b)H(a,b;x). (6) 

Here the integral (which may actually be a sum) covers 
all allowed a and 6, and the H (a, b; x) are the vector duals 
of the F(a, b; x). 

Wavelet transforms are superficially a generalization 
of Fourier transforms, for which the basis vectors are the 
exp(— ikx) and the dual vectors are the same function 
exp(ikx). In most cases, for wavelet transforms the F 
and the H differ markedly. Musical notation, now twelve 
centuries old, gives an early description of a continu- 
ous function using basic vectors-namely musical notes- 
that have local support, and that arise from translation 
and dilation. Equations I IKil assign to the function being 
transformed a continuous index x, but there is no math- 
ematical restriction that x be continuous. It is equally 
possible to apply wavelet transforms to a function gi that 
is labelled by a discrete index i. 

Wavelet and their transforms in general have a variety 
of important mathematical properties: 

(i) local support. The support of a function F is the re- 
gion over which it is nonzero. Most wavelets have only lo- 
cal support, and are non-zero over a limited range of x. In 
contrast, the Fourier transform function exp(— ikx) has 
non-local support, and is non-zero almost everywhere. 

(ii) local reconstruction. Changing a g{a, b) changes 
the reconstructed function (eqUJJ) over a limited region, 
namely the support of the associated H{a, b; x). 

(iii) On the other hand, most wavelet basis vectors are 
not antisymmetric or symmetric under x — > —x. Nor 
are they in general eigenfunctions of obvious differen- 
tial operators. Indeed, many wavelets only have a lim- 
ited number of derivatives. In some wavelet families, the 
F{a, 6; x) are not orthogonal; they are instead overcom- 
plete (though reconstruction is still possible). 

An extremely wide range of continuous and discrete 
wavelet transforms is known to exist. Much of the pub- 
lished literature on wavelet transforms examines one- 
dimensional transforms. To apply transforms to a volume 
one must either construct novel three-dimensional trans- 
forms, or one must construct a three-dimensional trans- 
form as the outer product of a series of one-dimensional 
transforms. We followed the latter approach. Space was 
divided into small volumes, a spherical harmonic den- 
sity was computed for each volume, and discrete wavelet 
transforms were then performed. 

We applied the discrete Haar^^ wavelets, which act 
by applying to adjoining pairs of points a low-pass filter 
C and a highpass filter D. The filters act on an adjoining 
pair of points £271-1 and x 2n as 

C„ = (x 2n -i + x 2n )/2 (7) 

and 

D n = {x 2n -i - x 2n )/2. (8) 
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The C filter produces a smoothed value, and the D filter 
acts as an edge detector. The transform is done itera- 
tively, by applying further pairs of Haar^^ filters to the 
output of the C filter from the previous iteration. Suc- 
cessive iterations are sensitive to features extending over 
larger spatial domains. To perform a Haar transform on 
a cubic array of points, one first applies the C and D fil- 
ters to pairs of points along one axis. The C and D filters 
are then applied to pairs of 1 x 2 blocks along the sec- 
ond axis, and then applied along the third axis to filtered 
2x2 plates along a third axis, so that, e.g., 

CCCm = Xiii+X 2 ll+£l21+£ll2+2;221+2;i22+£212+2;222 

(9) 

while 

DDDiu = 2111-2:211— £121— xn2+X22i+a;i22+a;2i2— £222- 

(10) 

Here Xijk is the value of x at the triple labelled point ijk. 
Because fluids have no directional orientation, the various 
CCD, CDC, and DCC filters should on thermal averag- 
ing give the same output. The wavelet decomposition is 
iterated by applying further cycles of C and D filters to 
the CCC output from the prior iteration, so that after 
each level of decomposition the CCC outcomes are aver- 
aged over a cubical volume with twice the linear extent. 
Physically, the CCC filter detects the uniformity in x 
across a region, the CCD filters are sensitive to surfaces 
bifurcating a region into zones with different values for 
x, the CDD filters are sensitive to lines where pairs of 
surfaces intersect, and the DDD filter identifies vertices 
where trios of surfaces come together. The filters other 
than CCC give signed output, so we computed averages 
such as {{CCD) 2 ). 



II. COMPUTATIONAL 

We studied a two-component A-B Lennard-Jones fluid. 
The potential energy was 



(ii) 



Here i and j label the species {A, B) of the interacting 
particles, so that Na and Nb are the number of A and 
B particles, while r is the distance between centers of 
mass. Parameters were taken from Glotzer, et al.fli). 
namely oaa = 1, Qbb = 5/6, and gab = 11/12. All 
particles had the same mass m and energy constant e. 
Natural units are used throughout, so the temperature 
T is units of fcs/e, m — 1 so force and acceleration are 
equal, and the time unit is a AA{m/ 'e) 1 ^ 2 ■ The potential 
was truncated: Uij(r) = for r > 2.5. The underlying 
molecular dynamics code is identical to that of our earlier 
paper 15] : Newton's equations of motion were integrated 
numerically using the Calvo and Sanz-Serna fourth order 
methodjla. Il7j. applying standard approaches to mini- 
mize calculation time. 



Our system is smaller (Na = 7813 as opposed to 
Na = 62, 500) than the system of our previous paper (l5|. 
We first did the calculations reported here and found that 
g{r) has a much larger range than sometimes reported. 
In order to determine g{r) accurately, the 15,625 parti- 
cle system of this paper was replaced with the 125,000 
particle system on which we have published 15]. The A:B 
mixture was almost exactly 1:1; Na and Nb differ by one 
atom. The density is the same as that of the previously- 
reported simulation: N/V = 1.30, where N is the total 
number of particles and V — L 3 is the system volume, 
with L — 22.96. Based on the range of known static 
correlations (at T = 0.56, g{r) -1/0 out to r > 10), 
the fluid studied here is weakly periodically confined, in 
the sense that the radial distribution shells around two 
particles can overlap twice, once in each of two opposite 
directions. 

We equilibrated our systems at temperatures 1.20, 
0.88, 0.73, 0.69, 0.66, 0.64, 0.59, and 0.56. The inter- 
mediate temperatures were chosen to match Ref. To 
equilibrate the system, we prepared an initial system at 
T = 1.20, equilibrated that system, and then cooled to 
new temperatures by rescaling the particle momentum by 
0.1% every fifty time steps, a time step being 0.01 in natu- 
ral units. On reaching the target temperature, the system 
was equilibrated for 10,000 time steps while monitoring 
the average kinetic energy. If the energy drifted during 
the 10,000 equilibration steps, the system was brought 
again to the target temperature and re-equilibrated until 
the average kinetic energy finally remained near its de- 
sired target value. Simulations at a given temperature 
were then run for at least 50,000 time steps. 

The generic radial distribution function g(r) (which 
gives the normalized (g(po) = 1) probability of finding 
any two atoms a distance r apart) and the specific radial 
distribution functions gij{r) (which give the probability 
of finding two particles of species i and j a distance r 
apart) were determined. The mean-square particle dis- 
placements, the (Q 2 ), and the Cjp{t) were then com- 
puted, taking the particles neighboring the particle of in- 
terest to be the particles in the first peak (region where 
g{r) > 1) of g{r), including I 6 (1,12). Static correla- 
tion functions {Q 2 } were evaluated using every particle 
as a center and repeating the calculation every fifty time 
steps. 

We examined a wavelet decomposition of several spher- 
ical harmonic component densities. In order to apply 
wavelet decompositions, we needed a rule for assigning 
wavelet densities to points ijk. We examined two rules, 
namely 
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Pijk ~ 
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Here the system is divided into cubelets having volumes 
Vijk, each containing atoms used as central particles 
for the evaluation of the Qg centered on them, and uq is 
the root mean square deviation of the distribution of the 
Qi. Equation 1121 gives the total deviation within Vijk of 
Qe, while eq EH gives the deviation as averaged over all 
the atoms within Vijk- Further work here is based on eq 

EH 

Finally, we examined a position-momentum time cor- 
relation function 



(YiP(t)) 



(Qio(t) Pz (t + T)) 



(14) 



(d QioW \ 2 )mw /2 

and a position-kinetic energy time correlation function 

(Q2o(t)K z (t + r)) 



(Y 2 K(t)) 



«|Q 20 (f) | 2 ><J3(t)»V3 



(15) 



Here p z is the ^-component of the momentum of the par- 
ticle of interest, while K z = pi /(2m), labels (t) and (t+r) 
refer to determinations on the same particle of interest 
at the indicated times, and the 9 — direction is taken 
to be the z-axis. For delay time t = 0, both (YiP(t)) 
and (Y^K^)) must vanish, because equal-time position- 
momentum correlations cannot exist in classical systems. 
However, in each correlation function the Ylm> Pi, and 
Ki were chosen so that they had exactly matching spa- 
tial tensor properties, so that the correlations in their 
products were not forced to zero by reflection or another 
spatial symmetry. For t ^ 0, we found that these corre- 
lations have a positive- or negative-going onset, respec- 
tively. 



III. RESULTS 

Figure^shows the dependence of (Qf) on I at temper- 
atures 0.56 and 0.88. The largest Qg are Q§, Q7, Q12, and 
Q 5 , in that order, followed by Qn and Qio- The signifi- 
cance of Q7 does not appear to be anticipated in the pub- 
lished literature, because no space-filling structure has 
7-fold (septahedral) symmetry, and the non-space-filling 
quasicrystalline symmetry that is usually discussed is the 
icosahedral symmetry, which emphasizes £ of 6, 10, and 
12. When the system is cooled from 0.88 to 0.55, Qq, 
and to a lesser extent Q7 and Q12, increase. All other 
Qg decrease with increasing temperature, the decreases 
in Q4, Q 5 , Q s , and Qg being the most pronounced. 

Figure shows the temperature dependence of the six 
largest (Qi). The (Qf) fall into three classes differing 
in their magnitudes and their temperature dependences. 
(Qg) is in the range 0.2-0.22; as the system is cooled 
from 0.88 to 0.56 it increases by almost 10%. (Q 2 ) and 
(Q12) are distinctly smaller than (Qg), namely they are 
rs 0.14 and increase modestly with decreasing temper- 
ature. (Qf), (Q 2 W ) and (Q 2 n ) are rj 0.1 - 0.12 at high 
temperature, and decrease by almost 10% as T is reduced 
from 1.2 to 0.56. (Q 2 n ) is sometimespj associated with 
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FIG. 1: Spherical invariants (Qf } as a function of I at temper- 
atures 0.55 (O) an d 0.88 (□). Note the marked amplitudes 
of the noncrystallographic £ = 5 and I = 7 invariants. 
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FIG. 2: Mean-square spherical invariants (Qf) as a function 
of temperature for £ of 5(0), 6(D), 7(a), 10(0), H(v), and 
12 (x). 



icosahedral ordering; we see it becomes smaller as tem- 
perature is reduced. 

We also determined the probability distributions, the 
likelihoods of determining particular values, for the (Qf). 
For each £, there are actually nine such functions, be- 
cause one can measure a (Qi ) averaged over all central 
particles, or only with an A or a B particle at the origin. 
Furthermore, in evaluating Qi around a given center, one 
could consider all particles in a given shell, or consider 
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only the A or the B particles in that shell. Limiting our- 
selves to an evaluation of Q 2 that includes all atoms in a 
given shell regardless of species, we found almost without 
exception that the probability distributions of the Q 2 are 
featureless, symmetric peaks whose widths at half height 
are typically half or 2/3 of the distribution's peak value. 
The exception is (Qf), as averaged over all central par- 
ticles, which has a prominent side shoulder, because the 
(Qg) for a central A and for a central B particle have 
very different average values. 

To further explore how the (Q 2 ) are affected by the lo- 
cal composition we determined the conditional averages 
{Qj)i,j,k, i.e., the value of (Q 2 ) for a given trio i,j,k. 
Here i is the identity of the central atom (A, B, All), 
while j and k are the number of A and B atoms, re- 
spectively, in the first shell; the calculation was made 
at T = 0.56. The number of particles instantaneously 
found in the first coordination shell fluctuates a great 
deal, from as few as four to as many (with an A par- 
ticle in the center) as fourteen. When plotted as func- 
tions of j and k, the (Qf)»,j,fc fell into three families: 
(i) Functions such as {Qf)A,j,k, (Qe)B,j,k, and (Qf 2 )s,i,)fc 
which decrease monotonically with increasing j + k, but 
are largely independent of j — k. (ii) Functions such as 
(Qs)B,j,k that have a minimum at intermediate j + k (for 
(Q\)B.j,k)i the minimum is at j + k = 8), that have max- 
ima at small and large j + k, and that are independent of 
j — k, and, finally, (iii) Functions such as {Q%)A,j,k that 
show a marked dependence on j — k, increasing at fixed 
j + k as the number of B atoms in the shell is increased. 

A qualitative explanation for the trends in (Q|)ij,fe 
with changing j and k may be suggested. As j + k is 
increased, the particles in the first coordination shell are 
obliged to pack more regularly rather than more ran- 
domly, at first making it more difficult for the instanta- 
neous Q 2 to assume large values. However, it is already 
quite difficult to pack 12 particles into the first coordi- 
nation shell of a B particle, as witness that for a B cen- 
ter atom the first shell practically never includes more 
than twelve atoms, so in order to obtain j + k = 12 the 
atoms in the first shell must be well-ordered, leading to 
a (Q1)bj,/c that increases at large j + k. Finally, the A 
particles are larger than the B particles, making it easier 
to fill the first coordination shell as the number of B par- 
ticles is increased, leading to a (Q|)A,j,fe that increases 
with increasing k — j. 

In addition to examining the static correlations of the 
Qi, we also examined their time autocorrelation func- 
tions. Preliminary studies found that all twelve G\ (t) 
relax on about the same time scale, so to reduce compu- 
tational demands we focused on the I having the largest 
magnitude, namely i — 5, 6, 7, 10, 11, and 12. Measured 
correlation functions at the highest and lowest temper- 
atures are found in Figure Figure^, at T = 1.20, 
represents a liquid that is warmer than the measured|l5| 
melting temperature T m w 1.1. At this temperature, the 
I = 10,11,12 fluctuations are nearly indistinguishable, 
and decay the most rapidly, the I = 6 and 7 correlations 
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FIG. 3: Spherically invariant autocorrelation functions 
C (2) {t) formed from the Flm(*) for fixed £ of 5(0), 6(»), 
7(B), 10(A), 11(#), and 12 (T) at temperatures (a) 1.20 and 
(b) 0.56. 



are the longest lived, and (t) has an intermediate life 
span. Each curve was fit to the sum of a short- and a 
long-lived relaxation, as described below. The fits are 
represented by the dashed lines, which at higher temper- 
atures are virtually indistinguishable from the data. The 
slow relaxations extend the relaxation curves to longer 
times, but are not apparent in the Figure as long-time 
shoulders. 

In FigureljJ), (T = 0.56, the lowest temperature stud- 
ied), the slow-mode shoulders are clearly apparent, while 



the Cg '(t) and C\ ' (t) functions are clearly distinct, 

(2) 

Cg (t) taking longer to decay. The amplitude of the slow 
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mode of C{ n (t) is markedly less than the amplitudes of 

(2) (2) 

the slow modes of C\ x (t) and C\ 2 (t). The accuracy with 

(21 

which the C\ (t) are described by the fitting process de- 
creases with decreasing T, so that in Figure^ the fitted 
curves (dashed lines) are no longer in near-perfect agree- 
ment with the measured data. 

In order to obtain a quantitative description of the 
relaxations, each relaxation function was fit separately 
to a form 

cf\t) = A s exp(-6 s t^) + A f exp(-8 f t f3 f). (16) 

Here subscripts s and / denote the slow and fast modes, 
Ai is an amplitude, 9i is a decay pseudorate, and Pi is a 
stretching exponent. A s and Af were taken to be inde- 
pendent variables, rather than being constrained to sum 
to unity. The Pi were constrained to lie in the interval 
[—1,1], but were free to vary within that interval (and 
in fact never entered the regime [-1, 0]), corresponding 
to processes that might be an exponential or a sum of 
exponential relaxations. Removing the constraint on the 
Pi sometimes gave very large values for Pi (1.4 to 2), with 
a very small increase in fit accuracy. When Pi = 1, the 
corresponding mode lifetime is 9~ 1 . In cases in which 
the relaxation was a stretched exponential in time, an 
average lifetime was assigned based on the exponential 
integral f^°dt exp(—9t) = 0~ x , namely 

poo 

(97 1 ) = / dtcxp(-9^) = 97 1/ ^T(1 + 1/ft), (17) 
Jo 

F representing the Gamma function. 

At low temperatures, best fits found pure exponential 
(P = 1) relaxations for both modes. At all temperatures, 
the shorter-lived mode was most often a pure exponcn- 

(2) (2) 

tial, though less so for Cg and C) than for the other 
relaxations. At higher temperatures, the slow mode is 
generally a stretched exponential, but as T declines pure- 

(2) 

exponential behavior is found. The slow modes of C{ , 

C^ii j and remain non-exponential to lower temper- 
atures than their lower-i? counterparts. 

Figures 0Ji and^J) show the lifetimes (9^ 1 ) of the fast 
and slow modes for each I at various temperatures. For 
i = 10,11,12, the fast mode lifetime increases gradu- 
ally with decreasing temperature, but not by more than 
30%. The slow mode lifetimes increase much more dra- 
matically than the fast mode lifetimes, by factors of 12 to 
20 over our temperature range. Both fast and slow mode 
lifetimes tend to converge with decreasing T, so that at 
T = 1.2 the lifetimes depend substantially on £, but at 
T = 0.56 the lifetimes of each mode are nearly indepen- 

(2) ' (2) 

dent of £. The fast-mode lifetimes of C§ , Cg , and 

(2) 

C) show a novel temperature behavior, namely a non- 
monotonic temperature dependence, with {9J 1 } climbing 
gradually to a peak at temperatures near 0.7, and then 
falling steeply as temperature is further reduced. 

Figures and^Jl give the amplitudes Af and A s for 
the fast and slow modes, respectively, against T. The 



two amplitudes were fit independently, and consistently 
give Af + A s « 1.1. The £ = 10, 11, and 12 relaxations 
are dominated by the fast mode, but for £ — 5,6,7 the 
slow mode has, at the least, nearly the same amplitude 
to the fast mode. The classes of temperature dependence 
seen above for the lifetimes are echoed by the amplitudes. 
For £ = 10, 11, 12, Af and A s are very nearly indepen- 
dent of £ and depend only weakly on T. The £ = 6 and 7 
modes show a non-monotonic temperature dependence, 
the reversal in slope being near T — 0.7. For £ — 5, the 
amplitudes are very nearly independent from tempera- 
ture. 

Figures^ and0f give the fast and slow mode lifetimes 
as a functions of £ for various temperatures. Figure 2t 
reinforces the prior description of mode behavior falling 
into three groups, namely £ = 5 (with a modest depen- 
dence of fast mode lifetime on T), £ = 6 and 7, with a 
wide range of (9J 1 ) with changing T, and £ = 10, 11, 12, 
whose behaviors are very similar to each other. Figure 0f 
presents the slow mode lifetime as a function of £. While 
at each £ there is a strong dependence of (9~ 1 ) on T, 
at fixed temperature the dependence of the slow mode 
lifetime on £ is much more modest. 

We now consider the wavelet decompositions of the 
spherical harmonic densities. We examined closely six 
values for £, three levels of decimation, and four filters, 
namely the mean-square outputs of the CCC, DCC, 
DDC, and DDD filters, all at nine temperatures. Be- 
cause the system is rotationally invariant, the DCC, 
CDC, and CCD filters are equivalent; we report the 
mean-square average of their outputs as DCC. Figures 
|5^ and|5j} show the spherical harmonic densities as func- 
tions of temperature at decimation levels of 1 and 3, i. 
e., the densities averaged over anSx5xS cube for S of 
2 and 8. From Figure decimation level 1, as T falls 
from 1.2 to 0.56 the smoothed densities for Q 5 and Qio 
decrease by nearly 30%, the smoothed densities for Qn 
and Q12 decrease slightly, and the smoothed densities for 
Qe and Qi increase, Q-j by more than Qq. On going to 
decimation level 2, which averages over 64 neighboring 
cubelets, the averaged Qe all decrease by nearly an order 
of magnitude, but the qualitative trends with tempera- 
ture remain the same. Between decimation level 2 and 
decimation level 3 (which averages over 512 neighbor- 
ing basic cubelets), leading to Figure 03, the smoothed 
Qi all fall by approximately another factor of five, and 
the temperature dependences of Q5 and Qg fade into the 
noise. Figure compares Q5 and Q7 at different dec- 
imations. The temperature dependence of Q5 decreases 
with increasing decimation, so that the first decimation 
of Q5 decreases markedly with decreasing T, but the sec- 
ond and third decimations depend rather weakly on T. 
In contrast, the temperature dependence of Q7 increases 
with increasing order of decimation, so that the second 
decimation increases more swiftly with decreasing T than 
does the first decimation, and the third decimation in- 
creases more swiftly with decreasing T than does the 
second decimation. 
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The numerical definitions of the DCC, DDC, and 
DDD filters are quite different from each other, but at 
each decimation the plots of their respective outputs, as 
functions of T for various t, are quite similar. The sim- 
ilarity arises because a change in temperature that in- 
creases the prominence of nominal vertices (as revealed 
by the DDD filter) simultaneously increases the promi- 
nence of edges and faces on a volume (the DDC and 
DCC filters), all by commensurate amounts. Figure 
shows the first and third ((DDD) 2 ) decimations as func- 
tions of T for various £. Just as the average CCC com- 
ponents of Q5 and Qio decrease as T falls (cf. FigJSJl, so 
also do the corresponding ((DDD) 2 ) components; simi- 
larly, just as the average CCC component of Q7 increases 
with falling T, so also do the corresponding ((DDD) 2 ) 
components. 

How are the changes in local statics and dynamics cor- 
related with other transport properties of the liquid? One 



approach to determining the local resistance of the fluid 
to flow is to examine the time-dependent self-diffusion 
coefficient D s (t) of an atom of fluid, as determined from 
the mean-square displacement of the component atoms 



((Ar(t)) 2 ) = 6tD a (t). 



(18) 



Here Ar(t) is the vector displacement of an atom during 
an interval t. D s (t) may in turn be used to compute a 
formal time-dependent viscosity rj(t) as 



D s (t) = 



k B T 
6m](t)R' 



(19) 



Here Boltzmann's constant ks is unity in natural units, 
T is the absolute temperature in the same units, and 
R is a particle radius. If the mean-square displacement 
of a Lennard- Jones fluid atom is used to compute D a 
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FIG. 4: Temperature dependence of the (a) fast and (b) 
slow mode lifetimes (0~ } and the (c) fast and (d) slow mode 
amplitudes A for fixed £ of 5(0), 6(D), 7(a), 10(0), H(v), 
and 12 (x). L-dependence of the (e) fast and (f) slow mode 
lifetimes at temperatures 0.56(0), 0.59(»), 0.62(D), 0.64(B), 
0.66(A), 0.69(0), 0.73(V), 0.88(+), and 1.20 (x). 



with a, b, c, and D s (oo) as fitting parameters. Up to 
constants, the long-time microscopic shear viscosity is 
~ T ' I D~ l , leading on identifying R = oaa to the nomi- 
nal microscopic viscosities seen in0). Over the observed 
temperature interval, the microviscosity changes by a fac- 
tor of 20, primarily at temperatures below 0.7. 



and thence rj, the inferred viscosity is properly a micro- 
viscosity whose relationship to the orthodoxly-measured 
macroscopic viscosity is complex. 

We evaluated ((Ar(i)) 2 ) for a randomly selected set of 
A particles. Representative determinations of ((Ar(i)) 2 ) 
against t appear in Figure^. At short times, ((Ar(i)) 2 ) 
has significant structure, but at long times the mean- 
square displacement tends to become linear in D s (oo)t. 
The measured displacements were fit to 

((Ar(f)) 2 ) = aexp(-rt c ) + D a (oo)t, (20) 



IV. DISCUSSION 

Here we have applied spherical harmonics and their 
invariants to study local orientational order in cold 
Lennard- Jones fluids. Time correlation functions and 
wavelet decompositions were employed to examine the 
temporal and spatial persistence of the orientation order. 
We began with a T = 1.20 fluid, a liquid slightly warmer 
than the apparent |l5j| crystalline melting point T m = 1.1, 
and finished with a deeply-cooled T = 0.56 vitrifying 
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fluid. Our particularly noteworthy findings include the 
importance of septahedral (£ = 7) order, the appearance 

(2) 

of slow modes in orientation relaxations C) (t) , and the 
relationship between the orientation fluctuation ampli- 
tudes, orientation fluctuations lifetimes, and the micro- 
viscosity. 

As seen in Figure |2 the largest (Qf) are those for i 
of 6, 7, and 12. (Qf) is next in importance, while (Q\^ 
is appreciably larger than (Q\q). Earlier literature on 
spherical harmonic components of the density largely re- 
stricted itself to even £, and thus did not identify the 
importance of (Qf) or (Q\x). 

The importance of (Qf) is readily understood; it re- 
flects a equatorial rings of six atoms surrounding the 
atom that is currently of interest. The large size of 
(Qj), which reflects the importance of seven- fold (septa- 
hedral)order, is unexpected and surprising. In our other 
work[lJ|, the most stable crystal structure that we found 
for this system was body-centered-cubic, which leads to 
amplitudes for (Qf) with £ even. No crystal structure has 
fundamental seven-fold orientation order. Indeed, seven- 
fold symmetry of any sort is uncommon in nature. A few 
complex organic compounds have distorted seven-fold 
symmetry [la. Il9| . as by forming heptamersj^. Com- 
puter simulations of heptamer rings have been madeplf. 
but no experimentally-observed quasicrystal has been re- 
ported to have seven-fold symmetry|22j. The identifiable 
importance of unidecahedral (£ = 11) ordering, as ob- 
served here, is also unexpected. 

The substantial sizes of the £ = 10 and £ = 5 compo- 
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FIG. 6: Temperature dependence of the (a) first and (b) 
third Haar wavelet ((DDD) 2 ) decimations of the spherical 
harmonic component densities Qe for I using symbols of Fig- 
ure 0] 



nents of (Qf) are less surprising. Proposals for the im- 
portance of quasicrystalline ordering in supercooled flu- 
ids |2^| have focused on icosahedral structures, which are 
not space-filling, but which are believed to be relatively 
stable for particles having Lennard- Jones interactions. A 
local fluctuation in the form of a body-centered icosahe- 
dron contributes to (Qi )- A pentagonal septode (a bi- 
furcated half-icosahedron) and fused pairs of icosahcdra 
also contribute to (Q\). 

The temperature dependence of the (Qf) speaks di- 
rectly to the potential importance of different aspects of 
orientational order in the transition toward the glass. As 
seen in Figure|3 as T falls, so do (Qf), (Q±o), and (Q n ); 
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in contrast, as T falls it is (Q\), (Qr), and (Q\ 2 ) that 
increase. Thus, to the extent that orientational ordering 
with £ = h or £ = \Q reflects the formation of icosahedra 
and pentagonal septodes, the presence of those structures 
decreases as the system becomes more viscous. Such a 
decrease in orientational ordering with decreasing T is 
contrary to any hypothesis that vitrification is associated 
with icosahedral ordering that is enhanced by reducing 
the temperature. The amplitudes that do increase are as- 
sociated with £ of 6 and 12, which arise with conventional 
close packing of spheres, and £ = 7, which describes an 
elsewise-unspecified but certainly not space-filling septa- 



hedral ordering. 

Instead of proposing that (Qfo) or (Q\ 2 ) is uniquely 
characteristic of icosahedral order, one might instead 
note that the £ — 11 component have properties very 
similar to that of the £ = 10 and 12 components: their 
slow and fast mode amplitudes are very nearly equal at all 
T, their fast mode amplitudes are approximately equal, 
and their {(CCC) 2 } and {{DDD) 2 ) wavelet components 
have about the same size and temperature dependences. 
It might then alternatively be said, based on our sample, 
that all large-£ components have similar static and dy- 
namic behaviors, and therefore the observed behavior of 
(Qi ) and (Q\ 2 ) does not appear to arise from tendencies 
for icosahedral clustering. 

Furthermore, there are considerable similarities be- 
tween the behaviors of the £ — 5, 6, 7 modes, similar- 
ities that make these modes systematically unlike the 
£ = 10,11,12 modes in their behavior, in particular: 
The slow and fast modes are similar in amplitude, the 
fast mode lifetimes depend non-monotonically on T, and 
for £ of 6 and 7 the mode amplitudes also depend non- 
monotonically on T, the slope reversals in all cases oc- 
curring near a T of 0.7. While the £ — 6, 10, 12 am- 
plitudes are associated with each other for icosahedral 
ordering, for the ordering actually found here it is the 
I = 7,6,5 amplitudes and separately the I = 10,11,12 
modes that appear to be associated with ordering in the 
fluid. Further study of higher-order spherical invariants 
might serve to clarify the nature of these orderings. 

Wavelet decompositions of the local spherical harmonic 
densities also speak to the relative importance of the vari- 
ous Qe . FigureJHt shows the first three decimations of the 
low-frequency ((CCC) 2 ) parts of the £ = 5 and £ = 7 den- 
sities at various temperatures. The amplitudes have been 
multiplied by constant factors, the same for all points of a 
given £ and decimation, so that fractional changes in the 
((CCC) 2 ) as T is changed are more readily compared. 
It is important to emphasize that each decimation sub- 
stantially reduces the ((CCC) 2 ), the reduction arising 
because with increasing decimation the ((CCC) 2 ) repre- 
sents an average of a fluctuation over a larger number of 
cubelets. In the Figure, this reduction has been masked 
by the multiplicative displacement. 

For £ = 5, the decline in the ((CCC) 2 ) with decreas- 
ing T is about the same for all three decimations, indi- 
cating that each decimation has about the same effect 
on ((CCC) 2 ) at all temperatures. On distance scales 
2 m aAA, m £ (1,3) probed by the decimations, the level 
of spatial coherence of (Q 2 ) as sampled by Haar low-pass 
wavelets is therefore independent of temperature. In con- 
trast, for £ = 7, the increase in ((CCC) 2 ) with decreasing 
T becomes more dramatic as the decimation level is in- 
creased. For the m = 1 decimation, ((CCC) 2 ) changes 
by 20% over the observed temperatures; for the m = 3 
decimation, the change in ((CCC) 2 ) is by nearly 2/3. Re- 
calling that each decimation greatly reduces ((CCC) 2 ), 
the actual observed effect is that, with decreasing T be- 
low approximately T = 0.7, decimations become decreas- 
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ingly effective in reducing ((CCC) 2 ), the mean-square 
fluctuation in Q-j as averaged over the region of sup- 
port of a Haar smoothing wavelet. Conversely, below 
T = 0.7, the range over which fluctuations in septahedral 
ordering are correlated begins to increase, out to at least 
the distances 2a aa to 8a aa over which Haar smoothing 
wavelets have support. To the extent that orientational 
ordering is significant for glass formation, it is only the 
spherical harmonic components that increase with de- 
creasing temperature that are plausibly contributing to 
vitrification, and these are the £ = 6 and I = 7 compo- 
nents but not the £ — 10 or 12 components. 

The ((DDD) 2 ) decompositions seen in Figure|S]mrther 
emphasize the possible importance of (Q\) and (Qf) in 
vitrification. For each decimation, {(DDD) 2 ) only in- 
creases with decreasing temperature for the £ = 6 and 
1 = 1 densities. For I = 5, 10, 11, and 12, at each decima- 
tion the high-frequency ((DDD) 2 ) components decrease 
or remain the same as T is reduced. 

The temporal persistence of fluctuations in the Qlm (t) 
is revealed by their time correlation functions, studied 
here from their spherical invariants C\ (t). The G\ (t) 
all have a bimodal structure, with a rapidly-decaying 
mode whose lifetime is nearly independent of T, and a 
slowly-decaying mode whose lifetime increases dramat- 
ically with decreasing T. The slow mode amplitude is 
quite weak at T = 1.2, is readily apparent at T = 0.88, 
and is much more prominent at lower temperatures. 

The fast mode has a relaxation time near 0.1 in natural 
units, relaxation times extending over no more than a fac- 
tor of 2, over the full range of T and £ that we examined. 
The relaxation time of the fast modes is very close to the 
periods of onset-time from initial zero to first maximum- 
of the correlation functions (Y±P(t)) and (Y2K z (t)). It 
might therefore be proposed that the fast mode corre- 
sponds to an initial period during which particle motions 
are highly correlated with aspects of initial particle posi- 
tions. 

The slow mode lifetimes increase from 0.7 ± 0.2 at 
T=1.20 to « 10 at T = 0.56. At higher temperatures, 
between £ = 6 and £ = 12, (9J 1 ) varies by by more than 
a factor of two, the larger-^ correlation functions hav- 
ing the shorter lifetimes. As temperature is reduced, the 
(O^ 1 ) converge toward a common £- independent value. 
The temperature at convergence is certainly not greatly 
different from the Mode Coupling Theory critical tem- 
perature estimated for this system by Sastry, et al.[2t|. 
We did not reach temperatures low enough to clarify if 
the relaxation times continue to keep a common value at 
temperatures colder than convergence, or if they again 
diverge from each other as T is further reduced. For 
£ = 10, 11, 12, the increase in log((6»7 1 }) with falling T is 
more rapid than linear in T, so that (9 J 1 ) for these modes 
might be diverging at some positive non-zero tempera- 
ture. In the direction of higher temperatures, extrapo- 
lation of the slow and fast mode lifetimes predicts that 
the lifetimes become equal near T ~ 2. This tempera- 
ture has had significance previously: In our previously- 



published[15( study of this system, a lower-temperature 
local structure formation was identified from the radial 
distribution function as appearing near T = 2 and be- 
coming more prominent as T is reduced. 

The ^-dependence of (9 J 1 ) suffices to prove that the 
relaxation of the C\ (t) cannot arise from the rotational 
diffusion of quasi-rigid evanescent structures, namely: If 
a rigid structure is performing rotational diffusion, the 
relaxation of its £ th spherical harmonic component will 
follow exp(— £(£+l)Dgt), where Dg is the rotational diffu- 
sion coefficient |24|. The relaxation time is inversely pro- 
portional to £(£ + 1), so for £ between 5 and 12 the time 
would change more than fivefold. The fast mode prop- 
erties are inconsistent with rotational diffusion: (9J 1 ) 
varies with £ by less than a factor of two; furthermore, 

(21 

the fast mode of C 12 !(t) is longer-lived than the same 

(2) 

mode of C\ x (t), which is in turn longer-lived than the 

( 2) 

fast mode of C[ (£), an order inverse to that expected 
from rotational diffusion. The slow mode properties are 
equally inconsistent with rigid-body rotational diffusion. 
Even at T — 1.2 the change in (9" 1 ) with £ is by less than 
a factor of three, is nearly independent of £ for £ in the 
range 10-12, and at low temperatures is very nearly inde- 
pendent of £ for all £ studied, entirely inconsistent with 
the rotational diffusion expectation that the rotational 
diffusion time should scale as £(£ + 1). 

Figure^ demonstrates that the apparent microviscos- 
ity of the liquid, as inferred from the particle self-diffusion 
coefficient, is correlated with the orientation amplitudes 
(Q|) and (Q 2 ). As indicated by the straight lines, for 
both amplitudes one sees 

T] = r lo cxp(a(Q 2 i )). (21) 

Here rj is rj at the perhaps- inaccessible (Q 2 ) = 0, and a 
is the slope seen in Figure^. Is this exponential correla- 
tion credible? There are other systems in which the vis- 
cosity increases more or less exponentially with the con- 
centration of the component responsible for the increase 
in viscosity. For example, in many polymer solutions 
over a wide polymer concentration range the viscosity of 
a polymer solution increases more or less exponentially in 
the polymer concentration!^. The correlation is there- 
fore not improbable on its face. From Figure |Ht>, the 
temperature dependence of the solution microviscosity is 
also not that different from the temperature dependences 
of the slow mode lifetimes. Over a shared range of de- 
creasing temperatures, rj and (9~ x ) each increase almost 
twenty- fold. 

There is a substantial literature on simulations of 
Lennard-Jones fluids, and a limited literature on spher- 
ical harmonic decompositions of particle densities. We 
first note several simulations that have found phenomena 
seemingly related to those studied here. We then turn to 
historical studies of spherical harmonic decompositions. 

For a system with the composition and potential en- 
ergy studied here, Sastry, et al.[2fi| found that the self 
intermediate scattering function F s a (q, t) of their A par- 
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FIG. 8: (a) Viscosity 77 from long-time self-diffusion coef- 
ficient D s (oo) plotted against (Ql) (O) an d (Qt) (•) an d 
(lines) fits to simple exponentials, (b) Viscosity (right axis, 
♦) and slow-mode relaxation times (left axis, points as per 
Figure as functions of temperature. 



tides gains a strong slow mode as T is taken from 0.88 
down to 0.73. The slow mode in F s a is a stretched expo- 
nential in time, whose stretching exponent (3 is « 0.9 at 
temperatures 0.88 and above; at temperature 0.73 and 
below (3 declines rapidly, reaching (3 k, 0.55 by temper- 
ature 0.59. The mode coupling critical temperature for 
the system was estimated to be T c = 0.592. Sastry, et al. 
also examined the distribution of particle displacements 
at fixed time and various temperatures, comparing distri- 



butions at the time at which the mean-square particle dis- 
placement was unity. The distribution of displacements 
was nearly Gaussian above T k, 0.63; at lower tempera- 
tures the distribution broadens and bifurcates. Between 
temperatures 0.88 and 0.59, the lifetime of the slow mode 
in F s A{q, t) increases by two orders of magnitude, consid- 
erably more than the increase that we found in (9 J 1 ) over 
the same range of temperatures. Nonetheless, the slow 
modes we find in C\ (t) are clearly accompanied by slow 
modes found in more familiar time correlation functions. 

A series of studies on Lennard- Jones fluids, e.g., by Do- 
nati, et al.j22}, have found that relatively mobile and rel- 
atively immobile particles separately tend to be found in 
clusters. The clusters of mobile particles are long and rib- 
bonlike; clusters of immobile particles are more compact. 
Donati, et al. present evidence that the liquid structure 
is much more ordered near an immobile particle than it is 
near a mobile particle, the additional ordering persisting 
out to a radius « aa- The distance over which immo- 
bile particles are correlated does not appear to increase 
much as temperature is reduced. On the other hand, the 
relative positions of highly mobile and highly immobile 
particles tend to be anticorrclatcd, the anticorrelation 
having a length scale that increases with decreasing T. 
From our wavelet decompositions, we find that the im- 
plicit range of some correlation functions, e.g., (<9§), does 
not appear to grow with falling temperature, just as Do- 
nati, et al., found that clusters of immobile particles do 
not appear to grow at T is reduced. On the other hand, 
we found from wavelet decompositions evidence that the 
range over which fluctuations in (Q7) are correlated ap- 
pears to increase as temperature is reduced, just as Do- 
nati, et al. found that range of anticorrclations between 
mobile and immobile particles grew as T was reduced. It 
will be interesting in future work to examine correlations 
between the details of local orientational order and the 
local mobility of particles of interest. 

Berthier [2^] determined a temperature-dependent 
range £ for spatial correlations of particle mobilities, 
showing that £ determines the distance scale at which a 
relaxation time r for particle motion decouples from the 
macroscopic viscosity 77, so that at lower temperatures rr] 
becomes temperature-dependent. Berthier proposes that 
his simulations show that the properties of supercooled 
liquids, and thus glass formation, are influenced by exten- 
sive spatial correlations, as seen previously in studies of 
critical phenomena. This proposal is said by Berthier[28) 
to be contrary to some models for glass formation. Fig- 
ure [S] shows a phenomenon similar to that discussed by 
Berthier, namely that the microviscosity (i.e., the self- 
diffusion coefficient) tracks accurately the degree of static 
orientational ordering for hexahedral and septahedral or- 
dering, and less precisely tracks the orientation fluctua- 

(2) 

tion correlation times from C\ (t). Our results differ 
from Berthier's in that our results speak directly to lo- 
calized ordering, involving the first coordination shell, 
rather than to the spatially extended ordering seen in his 
analysis. However, our wavelet decompositions indicate 
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that with declining temperature the septahedral order- 
ing, in addition to becoming more pronounced, is also 
increasing the range of its spatial correlations. Our anal- 
ysis does not show whether the relationship between the 
range increase we observe and the change in microscopic 
viscosity is causal or concomitant. 

The notion that local ordering in fluids may be non- 
crystallographic was proposed by Frank|llj. who noted 
that local icosahedral packing of Lennard- Jones particles 
was energetically preferred to hexagonal close-packed or 
face-centered cubic arrangements. Steinhardt, et al.0 
proposed to characterize local order in a fluid by making 
spherical harmonic decompositions of the local density of 
"bond" (near- neighbor vector) directions. They empha- 
sized the importance of spehrical invariants such as Qf. 
Icosahedral ordering contributes to invariants with £ of 
6 and 10, while cubic ordering has little invariant with 
£ = 10 but substantial invariant components with £ of 
4 and 8. Steinhardt, et al., report simulation results for 
even £ £ (2, 10), as well as for more complicated spheri- 
cal invariants, concluding that significant icosahedral and 
some cubic ordering is present. Tomida and Egami |5j ex- 
tended this approach by considering Q 2 t as computed for 
bonds centered on atoms, rather than spatial regions. 
Clusters centered on atoms provide a natural definition 
of local clustering. Their simulations used a potential 
suited for liquid iron, not the Lennard- Jones potential 
used here. They also found the probability distribution 
for the Q\. Clusters for which Q\ is large, as expected 
for icosahedra, were found by Tomida and Egami to have 
values for (Q2), (Q4), (Qg), an d (Q10) expected for icosa- 
hedra. However, the hypothesis of icosahedral ordering 
in fluids is not without controversy, as witness the careful 
analysis of Stillinger and LaViollete |29l l3Cl| . 

The results here do not entirely agree with the litera- 
ture, in part because we asked a somewhat different set 
of questions. In particular, we are not aware of prior 
measurements of (Qf) in this or similar systems, so prior 



studies would not readily have noticed the possible im- 
portance of septahedral ordering. Prior studies measured 
(Qio), finding that it was non-zero. Icosahedral, but not 
cubic, order leads to a large {Qio), leading to the con- 
clusion that icosahedral clustering is likely important in 
supercooled fluids. We measured not only (Qio) but also 
(Qil) and (Q12), finding that these three spherical har- 
monic components are very similar to each other in their 
magnitudes, relaxation times, and temperature depen- 
dences. Qn is not associated with icosahedral structures, 
so the observation that these three large-^ components 
have very similar behaviors suggests that the observed 
behavior of Qio and Q12 is some general large-^ effect, 
and not some consequence of icosahedral structure for- 
mation that would be expected to contribute to Qio but 
not to Qn- 

Kivelson, et al.|2^| have given a detailed model for 
glass formation in liquids, based on the formation of 
frustration-limited icosahedral clusters, starting at a clus- 
ter melting temperature significantly warmer than the 
melting temperature of the (unfrustrated) crystalline 
phase. In our previous paper, we proposed that we 
saw evidence for such a model, in the form of local, 
growth-limited clusters that appeared below a temper- 
ature T w 2 that is well above the crystalline melting 
temperature T rrn ~ 1.1 that we determined at our den- 
sity. One of us |31j has assembled evidence indicating 
that the so-called neutral polymer slow mode observed 
in light scattering spectra of some but not other nondi- 
lute solutions of neutral polymers provides experimen- 
tal evidence for Kivelson's model. In comparing results 
here with the proposals of Kivelson, et al., it is impor- 
tant to note that nothing in the Kivelson model relies 
on the growth-frustrated clusters being icosahedral. The 
non-space-filling septahedral (Qj) ordering for which ev- 
idence is found here is entirely consistent with the Kivel- 
son model's demands for a growth-limited spatial struc- 
ture. 
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